#!/usr/bin/python
import sys
import numpy as np

# parameters.
pileup_single_file = sys.argv[1]
pileup_double_file = sys.argv[2]
contig_size_file = sys.argv[3]
output_dir = sys.argv[4]

# Read in contig sizes.
print "reading contig sizes."
fin = open(contig_size_file, "rb")
lines = fin.readlines()
fin.close()

ctg_sizes = {}
for line in lines:
	tmp = line.strip().split()
	ctg_sizes[tmp[0]] = int(tmp[1])
del lines

# Build numpy arrays for each contig.
ctg_profiles = {}


# debug variable to detect most hit contigs.
debug = {}

# read in data.
print "reading single data."
fin = open(pileup_single_file, "rb")
lines = fin.readlines()
fin.close()

for line in lines:
	# tokenize.
	tmp = line.split()
	ctg = tmp[0]
	idx = int(tmp[1]) - 1
	cov = float(tmp[3])
	
	# check existance.
	if ctg not in ctg_profiles:
		tmp1 = np.zeros(ctg_sizes[ctg], dtype=np.float)
		tmp2 = np.zeros(ctg_sizes[ctg], dtype=np.float)
		ctg_profiles[ctg] = [tmp1, tmp2]
		
	# debug.
	if ctg not in debug:
		debug[ctg] = 0
	debug[ctg] += 1
		
	# add coverage.
	ctg_profiles[ctg][0][idx] = cov
del lines
	
# read in double data.
print "reading double data."
fin = open(pileup_double_file, "rb")
lines = fin.readlines()
fin.close()

for line in lines:
	# tokenize.
	tmp = line.split()
	ctg = tmp[0]
	idx = int(tmp[1]) - 1
	cov = float(tmp[3])

	# check existance.
	if ctg not in ctg_profiles:
		tmp1 = np.zeros(ctg_sizes[ctg], dtype=np.float)
		tmp2 = np.zeros(ctg_sizes[ctg], dtype=np.float)
		ctg_profiles[ctg] = [tmp1, tmp2]
	
	# debug.
	if ctg not in debug:
		debug[ctg] = 0
	debug[ctg] += 1
	
	# add coverage.
	ctg_profiles[ctg][1][idx] = cov
del lines

# scaled debug by 2 * size.
for ctg in debug:
	sz = 2 * ctg_sizes[ctg]
	perc = float(debug[ctg]) / float(sz)
	debug[ctg] = perc

# identify top 10.
tolist = sorted([(value,key) for (key,value) in debug.items()])
tolist.reverse()
toprint = []
if len(tolist) < 10:
	mx = len(tolist)
else:
	mx = 10
for i in range(mx):
	toprint.append(tolist[i][1])

print debug
print tolist

# print out profile.
print "printing output."
for ctg in toprint:
	# open output file.
	fout = open("%s/%s.csv" % (output_dir, ctg), "wb")
	for i in range(ctg_sizes[ctg]):
		if ctg_profiles[ctg][1][i] == 0.0:
			val = 0.0
		elif ctg_profiles[ctg][0][i] == 0.0:
			val = 0.0
		else:
			val = (ctg_profiles[ctg][1][i] / ctg_profiles[ctg][0][i])
			
		fout.write("%i\t%i\t%i\t%f\n" % (i, ctg_profiles[ctg][0][i], ctg_profiles[ctg][1][i], val))
	fout.close()
